rm(list = ls())
#setwd('~/Dropbox/Submitted Drafts/Compliance Blocking/PSRM Replication Files copy/Simulations')
sink("01_run_simulation.txt")
print(Sys.time())
start_time <- Sys.time()
print("Running Compliance Blocking simulations...")
options(dplyr.summarise.inform = FALSE)
library("quickblock")
library(tidyverse)
library(parallel)
library(estimatr)
print(sessionInfo())
source("Code/helper_functions.R")
#source("helper_functions.R")
#------------------------------------------------------------------------------------------------
n_iter = 5000
results_cov = c()
results_comp = c()
results_both = c()
set.seed(1)
for(n_sample in c(1000, 2000, 5000, 10000)){
  print(n_sample)
  results_cov = rbind(results_cov, mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = "covariates",
                         mc.cores = detectCores()) %>% bind_rows())
  results_comp = rbind(results_comp, mclapply(1:n_iter, run_iter,N= n_sample,  blocking_variables = "compliance",
                          mc.cores = detectCores()) %>% bind_rows())
  results_both = rbind(results_both, mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = c("covariates", "compliance"),
                          mc.cores = detectCores()) %>% bind_rows())
}
save(results_cov, results_comp, results_both, file="Simulation Results/sims_results_all.Rdata")
print("Completed Basic Simulations")
#------------------------------------------------------------------------------------------------
#ER Violations
 er_violation <- seq(0.1, 1, by =0.1)
 results_er = c()
 n_iter = 1000
 n_sample = 10000
 set.seed(1)
 for(b in er_violation){
   print(b)
   results_er = rbind(results_er,
   data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = "covariates", beta = b,
                       mc.cores = detectCores()) %>% bind_rows(), beta=b),
   data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = "compliance", beta = b,
                       mc.cores = detectCores()) %>% bind_rows(), beta=b),
   data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = c("covariates", "compliance"), beta = b,
                       mc.cores = detectCores()) %>% bind_rows(), beta=b))
 }
 save(results_er, file='Simulation Results/sims_er_violation.Rdata')
 print("Completed ER Violation Simulations")
#------------------------------------------------------------------------------------------------
#Compliance Probability:
 results_comp_prob = c()
 n_sample = 10000
 set.seed(1)
 for(i in c(0.1, 0.3, 0.5, 0.7, 0.9)){
  results_comp_prob = rbind(results_comp_prob,
                     data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = "covariates", compliance_rate = i,
                                         mc.cores = detectCores()) %>% bind_rows(), compliance_rate = i),
                     data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = "compliance", compliance_rate = i,
                                         mc.cores = detectCores()) %>% bind_rows(), compliance_rate = i),
                     data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = c("covariates", "compliance"), compliance_rate = i,
                                         mc.cores = detectCores()) %>% bind_rows(), compliance_rate = i))
 }
 save(results_comp_prob, file='Simulation Results/sims_comp_prob.Rdata')
 print("Completed Compliance Probability Simulations")
#------------------------------------------------------------------------------------------------
#PI Violations
pi_violation <- seq(0.1, 1, by =0.1)
results_pi = c()
n_sample = 10000
set.seed(1)
for(b in pi_violation){
  print(b)
  results_pi = rbind(results_pi,
                     data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = "covariates", gamma = b,
                                         mc.cores = detectCores()) %>% bind_rows(), beta=b),
                     data.frame(mclapply(1:n_iter, run_iter, N= n_sample, blocking_variables = "compliance", gamma = b,
                                         mc.cores = detectCores()) %>% bind_rows(), beta=b))
}
save(results_pi, file='Simulation Results/sims_pi_violation.Rdata')
print("Completed PI Violation Simulations")

print("Completed all simulations")
end_time = Sys.time()
print(end_time)
print(end_time - start_time)
sink()


